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Abstract 

The breakdown of finite element (FEM) computations for the flow of an Oldroyd-B 
fluid around a cylinder confined between parallel plates, at Weissenberg numbers 
Wi — 0(1), is shown to arise due to a coil-stretch transition experienced by polymer 
molecules traveling along the centerline in the wake of the cylinder. With increasing 
Wi, the coil- stretch transition leads to an unbounded growth in the stress maxi- 
mum in the cylinder wake. Finite element computations for a FENE-P fluid reveal 
that, although polymer molecules undergo a coil-stretch transition in the cylinder 
wake, the mean extension of the molecules saturates to a value close to the fully 
extended length, leading to bounded stresses with increasing Wi. The existence of 
a coil-stretch transition has been deduced by examining the behavior of ultra-dilute 
Oldroyd-B and FENE-P fluids. In this case, the solution along the centerline in the 
cylinder wake can be obtained exactly since the velocity field is uncoupled from the 
stress and conformation tensor fields. Estimation of the number of finite elements 
required to achieve convergence reveals the in-feasibility of obtaining solutions for 
the Oldroyd-B model for Wi > 1. 
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I. INTRODUCTION 



The behavior of polymeric hquids in complex flows is intimately linked to the distribu- 
tion of molecular conformations in the flow field. Macroscopic field variables such as the 
stress and velocity are strongly coupled to microscopic quantities such as the stretch and 
orientation of polymer molecules, and they influence and determine the magnitude of each 
other. Recent advances in computational rheology have led to the development of micro- 
macro methods that are capable of resolving information at various length and time scales. 
However, because of computational cost, most numerical simulations are still based on the 
purely macroscopic approach of continuum mechanics, where the conservation laws of mass 
and momentum are solved with a constitutive equation that relates t he stress t o the defor- 
mation history, without explicitly accounting for the microstructure Keunings 



2000j . The 



simplest constitutive equations capable of capturing some qualitative aspects of the vis- 
coelastic behavior of polymer solutions and melts are the Oldroyd-B and upper convected 
Maxwell models, respectively. These models have been widel y used in the investigation of 



Owens and Phillips 



2002 |. 



complex flows since the early days of computational rheology 

In spite of the apparent simplicity of the macroscopic equations, obtaining solutions at in- 
dustrially relevant values of the Weissenberg number Wi has proven to be extremely difficult 
in a range of flow geometries. Careful numerical studies over the past few decades suggest 
that the principal source of computational difficulties is the emergence of large stresses and 
stress gradients within narrow regions of the flow domain. Signiflcant efforts have been made 
to develop grid-based numerical techniques for resolving these stresses and their gradients. 
In spite of considerable progress, numerical solutions still breakdown at disappointingly low 
values of Wi ~ and it is still not clear whether this is because solutions do not exist at 

higher values of Wi, or w hether it is sim ply due to the in adequacy of current numerical tech- 



niques 



Keunings 



2000|. Very recently. 



Renardy 



2006| has shown analytically that in the 



special case of steady flows with an interior stagnation point, the mathematical structure of 
the upper convected Maxwell and Oldroyd-B models can be expected to lead to singularities 
in the viscoelastic stresses and their gradients with in creasing Wi. 



Nearly two decades ago. 



Rallison and Hinch 



1988( 1 argued in a seminal paper that in the 



case of the Oldroyd-B model, the inability to compute macroscopic flows at high Weissenberg 
numbers (the so called high Weissenberg number problem, or HWNP), has a physical origin 
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in a microscopic phenomenon. The Oldroyd-B model predicts an unbounded extensional vis- 
cosity in homogenous extensional flows at a critical value of Wi. The Oldroyd-B constitutive 
equation can be derived from kinetic theory by representing polymer molecules by Hookean 
dumbbells. The unphysical behavior in extensional flows is due to the inflnite extensibility 
of the Hookean spring used in the model. By considering the simple example of a stagnation 
point flow of an Oldroyd-B fluid, Rallison and Hinch showed that when the strain rate is 
supercritical, inflnite stresses can occur in the interior of a steady flow, brought about by 
the unbounded stretching of polymer molecules. Based on their analysis, they suggested 
the use of a constitutive equation that is derived from a microscopic model with a nonlinear 
spring force law (which would impose a flnite limit on a polymers extension), as an obvious 
re medy for the HWNP. 



Chilcott and Rallison 



1988| examined the benchmark complex flow problems of un- 



bounded flow around a cylinder and a sphere, using a dumbbell model with flnite exten- 
sibility, as a means of demonstrating the validity of this analysis. In order to understand 
the coupling between the polymer extension by flow, the stresses developed in the fluid, and 
the resultant flow fleld, they deliberately used the conformation tensor as the fundamental 
variable instead of the stress. The conformation tensor gives information on the distribution 
of polymer conformations within the flow fleld in an averaged sense. The use by Chilcott and 
Rallison of kinetic theory to develop their model enabled the derivation of a simple expres- 
sion relating the conformation tensor to the polymer contribution to the stress. By solving 
the equation for the conformation tensor along with the mass and momentum conservation 
laws, Chilcott and Rallison showed that even though there existed highly extended material 
close to the boundary and in the wake of the obstacle, there no longer was an upper limit 
to Wi in the range of values that was examined in their computations. Since the degree of 
molecular extension is directly related to the magnitude of stress, the Chilcott and Rallison 
procedure established a clear connection between high stresses and stress gradients in the 
flow domain with the conflgurational and spatial distribution of polymer conformations. In- 
deed, when simulations were carried out with the polymer length set to inflnity rather than 
a flnite value, the downstream structure was no longer resolvable, and the mean stretch of 
the polymers in the flow direction continued to grow with increasing Wi until the solution 
failed. 

In spite of this compelling demonstration of the physical origin of the HWNP, the up- 
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FIG. 1 Flow domain and boundary conditions for the flow around a cylinder confined between 
parallel plates. 

per convected Maxwell and Oldroyd-B models have continued to be used extensively in 
computational rheology. The reason for this might perhaps be attributed to the fact that 
even though stresses may be large, they are still bounded, and so far, there has been no 
conclusive demonstration that bounded solutions for the viscoelastic stress do not exist in 
complex flows at high values of Wi. On the contrary, by considering the steady flow of upper 
convected Maxwell and Oldroyd-B fluids around a cylinder confi ned between two parallel 
)lates (a schematic of the flow geometry is displayed in Fig. [1]), IWapperom and Renardy 
20051 ] have presented strong numerical evidence that suggests that solutions do exist for 
Wi > 1, and that current numerical techniques are not able to resolve them. 

The particular benchmark problem of flow around a cylinder between parallel plates was 
chosen by these authors since even though the geometry has no singularities, the maximum 
values of Wi for which converged solutions e xist are amongst the smallest of all benchmark 



flows 



Alves et al. 



2001 



Caola et al. 



2001 



Fan et al. 



199S 



Oscar et all 



2006 



Sun et al. 



1999| . The presence of upstream and downstream stagnation points leads to the development 
of steep stress boundary layers near the cylinder and in the wake of the cylinder, making the 
flow a stringent test of any numerical technique. Rather than solving the coupled problem 
simultaneously for the velocity and conformation tensor fields, for which the existence of 
solutions at high Wi is unknown, Wapperom and Renardy assumed a Newtonian-like velocity 
field and solved only for the conformation tensor field. The chosen velocity field has an 
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analytical representation, is very similar to the Newtonian velocity field in the same geometry 
(and the velocity field for an Oldroyd-B fiuid at Wi ~ 0{1)), and satisfies all the key 
requirements for the velocity field near the cylinder. The advantage of this approach is 
that the existence of a solution for an upper convected Maxwell model in such a velocity 



field, at all values of Wi is guaranteed Renardyl . |2000|, and hence failure of a numerical 



scheme can be attributed purely to numerics. By developing a Lagrangian technique which 
involves integrating the conformation tensor equation along streamlines using a predictor- 
corrector method, the authors were able to compute stresses up to arbitrarily large values 
of Wi (as high as 1024), and as a result, establish conclusively the existence of narrow 
regions with very high stresses near the cylinder, and in its wake. Further, they show 
that although the velocity field is known, one of the currently used numerical techniques, 
the backward-tracking Lagrangian technique, is unable to resolve the extremely thin stress 
boundary layers even for relatively low values of Wi. Since the fully coupled problem and 
the fixed fiow kinematics problem share the same basic dilemma of computing the stress 
field, Wapperom and Renardy argue that solutions also probably exist for Wi > 1 for the 
steady fiow of Oldroyd-B and upper convected Maxwell fluids around a cylinder confined 
between parallel plates, but current numerical techniques are not able to resolve them. 

The use of the conformation tensor as the fundamental quantity rather than the stress 
has become common in computational rheology, and the challenge of developing numerical 
methods capable of resolving steep stresses and stress gradients has been transformed to 
one of developing techniques capable of resolving rapidly varyirig confo rmation tensor fields. 
In an important recent breakthrough, iFattal and KupfermanI 2004| have shown that by 
changing the fundamental variable to the matrix logarithm of the conformation tensor, 
stable numerical solutions can be obtained at values of Wi significantly greater than ever 
obtained before. The success of their variable transformation protocol is predicated on their 
identification of the source of the HWNP as the inability of methods based on polynomial 
basis functions (such as finite element methods), to adequately represent the exponential 
profiles that emerge in conformational tensor fields in the vicinity of stagnation points and 
in regions of high defo rmation rate. 



Hulsen et al. 



2005l | have recently carried out a stringent test of the log conformation 



representation by examining the fiow of Oldroyd-B and Giesekus fiuids around a cylinder 
confined between parallel plates. (It is worth noting that unlike the Oldroyd-B models 
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prediction of unbounded extensional viscosity at a finite extension r a te, tlie extensional vis- 



cosity predicted by t 
fluids, Hulsen et al. 



le Gi esekus model is always finite 



Bird et al. 



1987aj). For both the 



2005| find that with the log conformation formulation, the solution 



remains numerically stable for values of Wi considerably greater than those obtained pre- 
viously with standard finite element (FEM) implementations. On the other hand, the two 
fluids differ significantly from each other with regard to the behavior of the convergence of 
solutions with mesh refinement. The lack of convergence with mesh refinement is usually 
seen most dramatically in the failure of different meshes to accurately predict the maximum 
that occurs in the normal polymeric stress component (Txx on the centerline in the wake of 
the cylinder. In the case of the Giesekus model, even at values of Weissenberg number as 
high as Wi = 100, mesh convergence is achieved in large parts of the flow domain, with the 
exception of localized regions near the stress maximum in the wake where convergence is 
not realized. For the Oldroyd-B model, however, the log conformation formulation fails to 
achieve mesh convergence in the entire wake region at roughly the same value {Wi > 0.6) 



as in previous studies. Further, the solution becomes unsteady at so me greater va . 



Hulsen et al. 



ue of 



2005 1 



Wi (depending on the mesh), finally breaking down at even higher Wi. 
speculate that this failure is probably due to the infinite extensibility of the Hookean dumb- 
bell model that underlies the Oldroyd-B model, and call for further investigations to see 
if this might lead to the non-existence of solutions beyond some value of Wi. Thus, after 
many years of attempting to resolve the HWNP by purely numerical means, it's origin still 
remains a mystery. 

In this paper, we conclusively establish the connection between the HWNP and the 
unphysical behavior of the Oldroyd-B model, in the benchmark problem of the steady flow 
around a cylinder confined between parallel plates. We show that when Wi ^ 1, polymer 
molecules flowing along the centerline in the wake of the cylinder undergo a coil-stretch 
transition, and that the location of the transition coincides with the maximum in the normal 
stress component axx on the centerline. With increasing Wi, the molecules stretch without 
bound, and this is accompan ied by a normal stress that increases without bound. (For an 



UCM fluid, 



Alves et al. 



200 1| have previously speculated that axx may develop a singularity 



at the position where it reaches a maximum value, based on the asymptotic behavior of highly 
accurate numerical results obtained by them using a finite volume method. However, they 
were unable to attribute a precise reason for the occurrence of a singularity). Computations 
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carried out here for a FENE-P fluid reveal that in this case also, polymer molecules undergo 
a coil-stretch transition which is located at the stress maximum. However, with increasing 
Wi, the mean extension of the molecules saturates to a value close to their fully extended 
length, enabling computations beyond the critical Weissenberg number. 

These dilute solution results have been obtained by drawing on insights gained from 
the solution of an ultra-dilute model, in which the velocity fleld is decoupled from the 
conformation tensor (an d stres s) fleld. This procedure is similar to the earlier work by 



Wapperom and Renardyl 2005| described above. However, rather than using an ad hoc 



velocity fleld, we have solved for the Newtonian velocity fleld using a full-fledged FEM 
simulation. Even for the ultra-dilute model, with a known velocity fleld, the solution of 
the equation for the conformation tensor (which is denoted here by M) for the Oldroyd- 
B model is found to breakdown for Wi ~ when a standard FEM method is used. 

However, an exact (numerical) solution for M, valid for arbitrary large values of Wi, is 
obtained along the centerline using two different techniques. In the flrst method, we exploit 
the fact that the equation for M in the Oldroyd-B and FENE-P models reduces at steady 
state to a system of ordinary differential equations (ODEs) along the centerline. In the 
second method, trajectories of the end-to-end vectors of an ensemble of dumbbells, flowing 
down the centerline in the wake of the cylinder, are calculated by carrying out Brownian 
dynamics simulations using the known velocity fleld. Averages carried out over the ensemble 
of trajectories lead to macroscopic predictions that are identical to the exact (numerical) 
results obtained by solving the macroscopic model ODEs discussed above, for arbitrary 
values of Wi. Comparison of the FEM results with the exact numerical results enables 
a careful examination of the reasons for the breakdown of the flnite element method. In 
particular, the occurrence of a coil-stretch transition at the location of the stress maximum 
in the wake is clearly demonstrated, and an estimation of the number of elements in the 
FEM method required to achieve convergence with increasing Wi reveals the in-feasibility 
of obtaining solutions for Wi > 1. 

The plan of this paper is as follows. In section [Tll we summarize the governing equa- 
tions, boundary conditions and computational method for the flow of dilute Oldroyd-B and 
FENE-P fluids around a cylinder conflned between parallel plates. We also elaborate on the 
connection between the macroscopic models, and the Kinetic theory models from which they 
are derived. In section UTTl we describe the means by which exact solutions to the ultra-dilute 
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models may be obtained, and examine the nature of the maximum in the Ma^^^ component 
of the conformation tensor. The resuhs of FEM computations and the exact numerical 
methods are first discussed for ultra-dilute solutions in section HVl followed by a discussion 
of the predictions of FEM computations for the dilute model. Section |V] summarizes the 
main conclusions of this work. 



II. DILUTE SOLUTIONS 
A. Basic equations 

As displayed in Fig. [1], the cylinder axis is in the z-direction perpendicular to the plane of 
flow. With the assumption of a plane of symmetry along the centreline {y = 0), computations 
are only carried out in half the domain. The cylinder, with radius a, is assumed to be placed 
exactly midway between the plates, which are separated from each other by a distance 2H. 
In common with other benchmark flow around a confined cylinder simulations, we set the 
blockage ratio H/a = 2. 

We normalize all macroscopic length scales with respect to a, velocities with respect to 
the mean inflow velocity far upstream <^f), macroscopic time scales with respect to 
and stresses and pressure with respect to ri(^v'^/a, where 1] = rjs + rj^ Q is the sum of the 
Newtonian solvent viscosity ?7s and the zero-shear rate polymer contribution to viscosity 
rjp^Q. Microscopic length and time scales are discussed subsequently. The two dimensionless 
numbers of relevance here are the Weissenberg number Wi = A(t')/a, in which A is a 
relaxation time, and the Reynolds number Re = pa (v') /rj, where p is the fluid density. 

The complete set of non-dimensional governing equations for a dilute polymer solution, 
described by the Oldroyd-B or FENE-P models, is 

V ■ V = (Mass balance) (1) 

Re V ■ Vv — Vp — V-Tg — V-cr = (Momentum balance) (2) 

^ + V ■ VM - K • M - M ■ = — ^ {/(tr M) M - 1} (Conformation tensor) (3) 

at Wi 

rs = 2/5D (Solvent stress) (4) 

cr = ^^^^^ {/(trM) M - 1} (Polymer stress) (5) 
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In these equations, k = (Vv)-^ is the transpose of the velocity gradient, D = | [i^'^ + i^) is 
the rate of deformation tensor, and the parameter /3 = {ris/r]) is the viscosity ratio. Here, we 
use (3 = 0.59, and set Re = 0, which are the values used in benchmarks for the Oldroyd-B 
model. The form of the function /(trM) depends on the microscopic model used to derive 
the equation for the evolution of the conformation tensor, and is consequently different in 
the Oldroyd-B and FENE-P models, as elaborated below. 



1. Oldroyd-B model 

As mentioned earlier, the Oldroyd-B constitutive equation can be derived from Kinetic 
theory, with the polymer molecule represent ed by a Hookean d umbbell model consisting of 



two beads connected together by a spring 



Bird et al. 



1987bl |. The spring obeys a linear 



spring force law F*^*^) = iJQ, where H is the spring constant, and Q is the connector vector 
between the beads. In this model, the conformation tensor M is defined by the expression, 

1 



M 



:qq> 



(6) 



where, ((.)) denotes an ensemble average, and {Q'^)^^ is the root mean square end-to-end 



vector at equilibrium. Note that, y (Q^)gq/3 = \/k-QT /H (with k-Q being the Boltzmanns 
constant and T the temperature), is the microscopic length scale. The connection between 
the microscopic and macroscopic models is established by deriving an evolution equation 
for the conformation tensor, and by relating the macroscopic polymeric stress to the confor- 
mation tensor. The expression of the conformation tensor is given by eqn ([3]) above, where 
in this case, the function / is identically equal to unity, z.e. /(trM) = 1. The polymer 



contribution to the stress is given by Kramers expression 

(T = nkjiT {M-I} 



Bird et al. 



1987b|, 



(7) 



where, n is the number density of polymer molecules in solution. With the relaxation time 
A, and the zero-shear rate polymer contribution to the viscosity r^p o, of the Oldroyd-B model 
related to microscopic model parameters by, 

c 



A 



AH 



Vp,o = nk^T A 



(8) 



where, C, is the bead friction coefficient, it is straightforward to see that the non- 
dimensionalization scheme used here leads to equation ([5]) for cr. 
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2. FENE-P model 



The FENE-P model corrects the faihng of the Hookean dumbbell model with its infinite 
extensibility, by using a spring force law that ensures that the magnitude of the end-to-end 
vector remains (on an average) below the fully extensible length of the spring Qq, 

H 



Q 



{i-{Q')/Ql) 

For the FENE-P model, the conformation tensor is also defined by eqn 
mean square end-to-end vector at equilibrium (Q^)^^ given by, 

3QHkBT/H) 



eq 



Ql + 3 {ksT/H) 
Defining the finite extensibility parameter b by, 

, Ql 



(9) 

but with the 
(10) 



(11) 



eq 



the function /(trM) for the FENE- P model in the confo r matio n tensor evolution equa- 



tion ([3]) can be shown to be given by 



Pasquali and Scrivenl . 



2002j . 



/(trM) 



(12) 



b-ti M/3 

Note that the definition of b us ed in eqn (ITTll is d ifferent from the definition of the finite 



extensibility parameter given in 



Bird et al. 



1987b|, which is widely used in the literature. 



In the Bird et al. definition, the microscopic length scale used to non-dimensionalize Qq is 
the Hookean dumbbell length scale \Jk^T/ H. The polymeric stress for a FENE-P fiuid is 
given by, 

<T = nA;BT {/(trM)M-I} (13) 

With the relaxation time A, and the zero-shear rate polymer contribution to the viscosity 
?7p_o related to microscopic model parameters by, 

'b-l 



A 



AH 



Vpfi = nkjiT A 



(14) 



non-dimensionalization of eqn (fT3l) leads to eqn 
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TABLE I Meshes used in the finite element simulations for computing viscoelastic flow around a 
cylinder. 



Mesh 


Ml 


M2 


M3 


M4 


M5 


Number of Elements 


2311 


5040 


8512 


16425 


31500 




FIG. 2 Mesh M2 used in the finite element simulations. 



B. Boundary conditions and computational method 



m 



by 



The set of governing equations ([I])-© are solved with the boundary conditions shown 



le location of the inflow and outflow boundaries coincides with that chosen 



Sun et al. 



1999| . who showed that the flow is insensitive to further displacement of 
the open boundaries in the range of Weissenberg numbers examined. A no-slip boundary 
condition is imposed on the cylinder surface and the channel walls. Fully developed flow 
is assumed at both inflow and outflow boundaries, with the velocity prescribed at both 
boundaries. Here, (v'^ = 1 for the prescribed velocity field. The boundary conditions on 
the conformation tensor are imposed only at the inflow boundary. At the symmetry line, 



tn : fr. 



and Vy = is imposed, where, t and n are the unit vectors tangential and 



normal to the symmetry line, respectively. 

The governin g equations are discretized by using the DEVSS-TG/SUPG mixed finite 



element method Pasquali and Scriven 



2002| . The DEVSS-TG formulation involves the in- 



troduction of an additional variable, the velocity gradient k. Continuous biquadratic basis 
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functions are used to represent velocity, linear discontinuous basis functions to represent 
pressure and continuous bilinear basis functions are used for the interpolated velocity gradi- 
ent and conformation tensor. The DEVSS-TG/SUPG spatial discretization results in a large 
set of coupled non-linear algebraic equations, which are solved by Newton's method with an - 



alytical Jacobian and first order arc-length continuation in Wi [Pasquali and Scrivenl . 12002 1. 
Five different meshes are used for the FEM calculations. Details of the different meshes used 
in this work are given in Table [H and the mesh M2 is displayed in Fig. [2j The important 
distinction among the five meshes is the density of elements on the cylinder surface and in 
the wake of the cylinder. 

III. ULTRA-DILUTE SOLUTIONS 



The phrase "ultra-dilute" is used to describe a situation where, even though polymer 
molecules are prese nt, they have a negligible effect on the velocity field. As demonstrated in 



the earlier work by IWapperom and Renardyi 2005| , the solution of the ultra-dilute problem 



provides important insights into the structure of the solution to the dilute problem, where 
the polymer stress and velocity fields are fully coupled. As will be clear from the results 
and discussion presented subsequently, the solution of the ultra-dilute case lies at the heart 
of the analysis carried out in this work. 

Since the velocity field for an ultra-dilute solution is determined completely by the solvent 
stress, it is identical to the velocity field for a Newtonian fluid. The ultra-dilute conformation 
tensor and velocity fields are simply obtained by solving the governing equations with the 
parameter /3 set equal to unity. As can be seen from eqn ([5]), this implies cr = 0, leading 
to eqns and ([2]) being identical to the mass and momentum balances for a Newtonian 
fluid. The same FEM formulation described above for the solution of the dilute case, can 
consequently, also be used to obtain the ultra-dilute conformation tensor and Newtonian 
velocity fields when (3 = 1. 

The conformation tensor field, for an ultra-dilute solution, corresponds to the average 
configurations of polymer molecules in a pre-determined velocity field. Note that the choice 
of value for Wi has no influence on the velocity field (which is determined once and for all 
for the specified geometry), but can significantly effect the conformation tensor field, as will 
be seen subsequently. 
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Clearly, the polymer contribution to the stress tensor cannot be obtained from a solution 
of the macroscopic equations, where it is assumed to be zero. In this case, we resort to 
reporting the dimensionless stress predicted by the microscopic models, cr* = cr/nk^T = 
(A/r/p o) cr, given by eqns ([7]) and f[T^ . in the Oldroyd-B and FENE-P models, respectively. 
This is similar to the evaluation of the polymer contribution to the stress in a homogenous 
simple shear or extensional flow, where the velocity field is prescribed a priori. 

As mentioned earlier, the FEM formulation used here fails to give mesh-converged results 
for the conformation tensor field for an ultra-dilute solution beyond a threshold value of the 
Weissenberg number. The reasons for this failure are analyzed subsequently. Crucial for 
this analysis, however, is the possibility of obtaining an exact (numerical) solution to the 
conformation tensor equation along the centerline in the wake of the cylinder. Two means 
of obtaining such a solution are discussed below. 

A. Exact (numerical) solution along the centerline in the cylinder wake 

1. System of ODEs 

By the requirements of symmetry along the centerline, the velocity field must have the 
form, Vx = g{x), Vy = 0, and the components of the velocity gradient tensor must satisfy, 
i^xy = i^yx = 0. Incompressibility requires k^x = ~i^yy The flow along the centerline is 
consequently planar extensional in character. 

Substituting these results into the evolution equation for the conformation ten- 
sor (eqn ([3])) leads, at steady state, to a system of ODEs in the independent variable x, 
for the components of the conformation tensor M. In the case of the Oldroyd-B model, the 
equations for each of the components are decoupled from each other. For the purposes of 
analysis in the present work, we are only interested in the equation for the M^x component, 
which can be shown to be. 



Since the velocity field is known a priori, equation ( |T5i) is a first order liner ODE for 
Mxx that can be solved straight-forwardly once a boundary condition is prescribed. We 
can expect that at the downstream stagnation point at x = 1, the Hookean dumbbells that 
represent the polymer molecules are in their equilibrium configurations, and consequently 




(15) 
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Mxx = 1- However, as can be seen from the form of eqn (fTSl) . we cannot use this boundary 
condition since f = at the stagnation point. We can overcome this difficulty by exploiting 
the fact that we can guess the form of the velocity field asymptotically close to the stagnation 
point. For an unbounded flow of a Newtonian fluid near a stagnation point on the surface of 
a two-dimensional body, the assumption of a quadratic velocity fleld, asymptotically close 
to the stagnation point, is exactly v alid for Stokes flo w, and is consistent with the accepted 



numerical solution for Hiemenz flow 



Pozrikidid. 



fleld postulated by IWapperom and Renardv 



also of the form Vx = k {x — 1) (see also [Renardv 



2005 



19971]. The approximate Newtonian velocity 



'or th e flow around a conflned cylinder is 



20001), with = 4, in the limit 



(from above). We flnd that the assumption of a quadratic velocity fleld, with a value of 
k = 4.178, leads to an excellent flt of the Newtonian velocity fleld close to the downstream 
stagnation point obtained by the FEM solution. For a quadratic velocity fleld, eqn f|T5|) 
admits an analytical solution for M^x, 

M^^ = 1 + 2a {x-l) + Sa^ {x - if + {x - if + 1.5a^ {x - 1)^ (16) 

where, a = 2kX. As a result, we use the analytical value of Mxx at x = 1.01 as the boundary 
condition to integrate eqn (fT5|) . with a Runge-Kutta 4th order method. The functions Vx{x) 
and obtained by interpolation from the Newtonian FEM solution, at each of the 

values of x where they are required for the purpose of integration. 

For the FENE-P model, the diagonal components are not decoupled from each other, and 
consequently, evaluation of the M^^^ component requires a solution of a system of ODEs for 
all the diagonal components, as can be seen from the equations below, 

2\Kxx{x) > Mxx + 



dx y^Vx{x) [6 — trM/3 j Af^(x) 

2Ak,,w}m„ + ^ (17) 



dUyy _ 1 r 6-1 , A 1 



dx \vx{x) [6 — trM/3 

dM,, 1 r 6-1 



dx Xvx{x) 1^6 — trM/3j Xvx{x) 

It is difficult to solve these equations analytically even with the assumption of a quadratic 
velocity field close to the stagnation point. Since i^xx{,x) ^ 1 at x = 1.01, we use the equilib- 
rium initial conditions, Vixx = 1, = 1, and Mzz = 1, as the boundary conditions, though 
these values are strictly correct only at the stagnation point x = 1. It turns out, however, 
that the conformational tensor fields along the centerline, downstream of the stagnation 
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point, are insensitive to a variation by a few percent, in the boundary values chosen for the 
diagonal components of M close to the stagnation point. (For instance, identical results are 
obtained if the boundary conditions for the Oldroyd-B model at x = 1.01, are used instead). 



2. Brownian dynamics simulations 

An alternative means of obtaining an exact solution along the centerline in the cylinder 
wake is to exploit the connection between the macroscopic and microscopic models. If 
we imagine an ensemble of dumbbells at any position x on the centerline, subject to the 
local velocity gradient, we expect that at steady state, the ensemble average (Q^Q^) = M, 



where, Qt = Q/ ^i{Q'^)jS). We could also imagine a packet of fluid with an ensemble 
of dumbbells, starting close to the stagnation point and traveling down the centerline with 
velocity Vx{x), experiencing the local velocity gradient at each position x. In this case, the 
variation with time of (Q^Q^) would be equivalent to the variation of M with x in the 
macroscopic models. The ensemble average (Q^Q^) can be obtained by integrating the 
stochastic differential equation (SDE), 

dQt = Mt) -Q^-^J {{Q'')) Qn dt + ^^ dW^ (18) 

which governs the stochastic dyna mics of Hookea n or FENE-P dumbbells subject to the 



Ottinger 



19961]. Here, is a non-dimensional Wiener 



time varying velocity gradient /t(t) 
process, and 

/ r, N 1 for Hookean dumbbells, 

f{{Q')) = { , (19) 

I (6 - l)/{b - (Q^ >/3) for FENE-P dumbbells 

is the same function for the FENE-P model as in eqn f[T^ . with (Q^^) taking the place of 
trM. 

The only remaining issue is to obtain K,xx{t) for a packet of fluid traveling down the 
centerline in the wake of the cylinder. This can be done in a straightforward way since we 
know Vx{x) and k,xx{x) for a Newtonian fluid. Since dx/dt = Vx{x), the integral, 

t = [ dx' = hix) (20) 

Ji+s Vx[x') 

gives us t as a function of x for a material particle. Clearly, K,xx{t) = i^xx{h~^{t)). The 
fluid packet would have to start its journey slightly downstream of the stagnation point 
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1 5 10 15 20 25 30 

t 

FIG. 3 (Color online) Velocity {vx) and velocity gradient (k^x) for an ultra-dilute solution (or 
a Newtonian fluid) along the center line in the wake of the cylinder, computed using the FEM 
formulation with the M5 mesh. The dashed green line is the position x, as a function of time, of 
a material particle traveling downstream starting close to the stagnation point. The dashed black 
line is the time dependent velocity gradient Kxx{t) used for carrying out BDS of the ultra-dilute 
models. The dot-dashed line is the velocity profile for a dilute Oldroyd-B model at Wi = 0.6. 

(represented by 1 + 5 in the lower limit of the integral in eqn (!20l) ). since otherwise it would 
remain indefinitely at the stagnation point. Time is therefore initialized at x = 1 + 5. 

The velocity field Vx{x) and the velocity gradient k,xx{x) along the centerline, computed 
using the FEM formulation with the M5 mesh, are displayed in Fig. [3l Note the quadratic 
nature of the velocity close to the stagnation point. The position x of a fiuid packet as a 
function of time t, calculated using eqn (120|) . is displayed as the dashed green line. Consistent 
with the boundary conditions for the ODE for Mxx above, we use 6 = 0.01. Clearly, a 
material particle spends a significant fraction of time close to the stagnation point before 
rapidly accelerating away from it. The time dependent velocity gradient Kxx{t), necessary 
for the integration of the stochastic differential equation f|T8l) . calculated as described above, 
is shown as the dashed black line. 

The configurational distribution function for the Hookean dumbbell model is Gaussian 
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both at equilibrium, and in the presence of a homogenous flow fleld [Bird et all Il987b |. 
Since a Gaussian distribution is completely determined by its second moments, the initial 
distribution function at x = 1.01 can be calculated for the known velocity gradient k^j.^, 
using the analytical solution for M^^, in eqn flTBl) (and a similar expression for Myy. Note 
= 1). For Hookean dumbbells, the SDE (ITS]) is integrated forward in time her e, using a 
secon d order predictor-corrector Brownian dynamics simulation (BDS) algorithm Ottingerl . 



1996|, with an initial ensemble of connector vectors distributed according to the Gaussian 



distribution at x = 1.01, subjected to the time dependent velocity gradient Kxx{t) for t > 0. 

The distribution function is also Gaussian in the case of the FENE-P model, and as a 
result, the initial distribution function at x = 1.01 can in principle be determined using 
the second moments obtained by solving the set of governing equations f[T7|) . However, we 
adopt the simpler procedure of using a Gaussian distributed initial ensemble with equilibrium 
second moments, M^^ = 1, Myy = 1, and M^z = 1, at x = 1.01, since the solution of the 
SDE downstream of the stagnation point is found to be insensitive to the choice of the 
initial distribution. (For instance, identical results are obtained if the initial distribution 
of Hookean dumbbells at x = 1.01, is used instead). The BDS algorithm for FENE-P 
dumbbells is identical to the one used to integrate the SDE for Hookean dumbbells, with 
the additional feature of having to evaluate (Q^^) at every time step. 

The exact (numerical) results along the centerline in the cylinder wake, for an ultra-dilute 
solution, calculated by solving the ODEs and the SDE above, are compared with the FEM 
solution in section [IV] below. Before doing so, however, important insights can be obtained 
by considering the nature of the maximum in the M^x component in the wake of the cylinder. 



B. The maximum in the wake 



For an ensemble of dumbbells which start near the stagnation point (where the velocity 
gradient gligibly small), and then travel downstream to the region of fully developed 

flow (where k^x = 0), the M^^, component of the conformation tensor initially has a value 
close to unity, and ultimately returns to a value of unity. Since Kxx > at intermediate 
values of x, it is clear that M^^^; must attain a maximum at some point x = x* along the 
centerline in the wake of the cylinder. Indeed all computations of flow around a conflned 
cylinder report the occurrence of a maximum, and as mentioned earlier, failure to attain 
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FIG. 4 (Color online) The dimensionless polymer contribution to stress, cr*^ = {X/r]pfi)axx, along 
the cylinder wall and along the centerline in the cylinder wake, at Wi = 1.0, for an ultra-dilute 
solution of an Oldroyd-B fluid. The computations are carried out using the FEM formulation 
discussed in section Ill.BI The inset shows the lack of mesh convergence in the wake. 

mesh convergence is usually observed most noticeably at the maximum. 
At the maximum (where dM^^/dx = 0), eqn 0151) implies that, 

XX x^^x * 

Clearly, if Xkxx = 0.5 at x = x*, the maximum value of M^x will be unbounded. One can see 
from the dotted curve in Fig. [3] that for values of A = 0{1), there are many points along the 
centerline in the cylinder wake where Xkxx can be greater than 0.5. However, as manifest 
from eqn fl?I]) . the real issue is whether Xkxx = 0.5 at x = x*. This question is examined in 
the section below, using the different solutions methods discussed above. 

IV. RESULTS AND DISCUSSION 
A. Ultra-dilute solutions 

The failure to attain mesh convergence in the cylinder wake, for an ultra-dilute polymer 
solution at Wi = 1, using the FEM formulation discussed in section III.Bl is displayed in 
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terms of the non-dimensional polymer contribution to the stress in Fig. HI The profile for 
the stress, with two maxima, one on the cylinder wall, and a second in the wake is typical 
for viscoelastic flow around a confined cylinder. As Wi increases further, the maximum in 
the wake grows significantly, and becomes the more dominant of the two maxima. The 
lack of mesh convergence, which is already apparent at Wi = 1 in Fig. HJ becomes much 
more pronounced. These predictions are completely in accord with what has been observed 
previously for dilute solutions. They have only been reproduced here to demonstrate the 
existence of a similar mesh convergence problem even in the simpler case of an ultra-dilute 
solution. 

The advantage of considering ultra-dilute solutions is that exact (numerical) solutions 
along the centerline in the cylinder wake can be obtained as outlined in section IIII.AI Fig. [S] 
compares the predictions of FEM computations with the exact solutions obtained by solving 
the ODEs, eqns (fT5ll and (fTTIl . and by integrating (using BDS) the SDE, eqn (fT8|) . for the 
Oldroyd-B and FENE-P models. The excellent agreement of the FEM results with the 
exact solution (except for the Oldroyd-B model at Wi = 1.3) shows that the FEM results 
are accurate at the Weissenberg numbers that have been displayed. The departure of FEM 
predictions, at Wi = 1.3, from the exact results for the Oldroyd-B model (Fig. [5] (b)), is 
clear proof of the breakdown of FEM computations for Wi > 1. Since the exact solution 
is known for any value of Wi, the error in the FEM computation can be estimated. Before 
we discuss the error, however, we first consider the more pressing issue of the value of the 
non-dimensional strain rate Xkxx a.t x = x*, the location of the maximum in M^x in the 
cylinder wake. 

Fig. [6] (a) displays Xkxx\x=x* as a function of Wi, for an ultra-dilute Oldroyd-B fluid, ob- 
tained with the three different solution techniques. As was discussed earlier in section IIII.B[ 
the maximum in M^^^^ becomes unbounded as Xk,xx\x=x* 0.5 (see eqn (12T!) ). Interestingly, 
Xkxx\x=x* first approaches 0.5 at Wi ~ 1, where computational difficulties with the FEM 
method are first encountered. While the ODE solution and BDS can be continued to higher 
Wi, FEM computations (indicated by the crosses) are no longer accurate beyond Wi = 1, 
breaking down completely by Wi = 1.55. This is related, as will be discussed in greater 
detail shortly, to the inability of the FEM solution to resolve the very large stresses that 
arise as Xkxx\x=x* approaches 0.5. Fig. [6] (b) shows that the approach of Xkxx\x=x* to the 
critical value is approximately linear in Weissenberg number for small values of Wi, becom- 
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(c) Wi = 1.5 (d) Wi = 2.0 



FIG. 5 (Color online) The M^.^. component of the conformation tensor, computed using FEM, 
ODE and BDS, for an ultra-dilute solution of Oldroyd-B and FENE-P (6 = 55) fluids. The EDS 
results are obtained by averaging over 10^ dumbbell trajectories. The FEM results are computed 
using the M5 mesh. For the Oldroyd-B model, at Wi = 1.3, the FEM results at the maximum are 
approximately 6% different from the ODE or BDS values. 

ing a power-law (Wi~'^'^^), for Wi > 0.4. Thus, though Xkxx\x=x* — > 0.5 asymptotically, it 
never equals or exceeds the critical value. This implies that M^x (and, consequently, a*^) 
will increase without bound as Wi increases, but will never become singular. 

In the case of the FENE-P model, Xk,xx\x=x* approaches and exceeds 0.5 with increasing 
Wi, as shown in Fig. [71 for a range of values of the finite-extensibility parameter b. The 
value 0.5 is not significant for the FENE-P model since, as can be seen from eqn (fTTI) . 
Mxx\x=x* — > (& — trM/3)/(trM/3 — 1) as Xk,xx\x=x* 0.5. Hence, Mxx\x=x* remains 
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Wi log Wi 

(a) (b) 

FIG. 6 (Color online) (a) Dependence of the non-dimensional strain rate Xkxx on Wi, at the location 
X = X* the maximum in yixx in the cylinder wake, for an ultra-dilute Oldroyd-B fluid. Inset 
shows \k,xx\x=x* approaching 0.5, computed from the ODE solution, for Wi > 1.5. (b) Xkxx\x=x* 
approaches the critical value 0.5 as a power-law with increasing Wi. The FEM results are for the 
M5 mesh and BDS results are obtained by averaging over 10^ individual Brownian trajectories of 
Hookean dumbbells. 

finite as long as b is finite. Before we compare the predictions of Mxx\x=x* and a*^\x=x* by 
the Oldroyd-B and FENE-P models as a function of Wi, it is interesting to consider the 
dependence of the location of the maximum, x = x*, and Kxx\x=x*, on Wi. 

Figure [S] (a) shows that for the Oldroyd-B model, the location x* of the maximum in Mxx 
continuously moves downstream away from the stagnation point in the cylinder wake, with 
increasing Wi. Simultaneously, Fig. [8] (b) indicates that the value of Kxx\x=x* continuously 
decreases. This is consistent with the behavior of function of x displayed in Fig. [31 

Since the velocity field is determined a priori for an ultra-dilute solution, Wi is varied 
in the computations by varying A. With increasing Wi, the product \kxx\x=x* (with A 
increasing, and Kxx\x=x* decreasing) tends to 0.5 in the manner depicted above in Fig. [n](a). 
The continued use of the Oldroyd-B fluid in computational rheology, in spite of its known 
shortcoming in extensional flow, has sometimes been justifled by the argument that real flows 
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FIG. 7 (Color online) Dependence of the non-dimensional strain rate Ak^xj at the location x = x* 
of the maximum in M^..^ in the cylinder wake, on Wi, for a FENE-P fluid. 

adapt to avoid an infinite stress. Curiously, the results in Figs. [6] and M appear to suggest 
that the flow does regulate the stress, and avoid a singularity. This does not, however, 
prevent the failure of the FEM computations for the Oldroyd-B fluid! 

For a FENE-P fluid, the location of the maximum in M^x and the value of Kxx at x = 
X*, display intriguing behavior with increasing Wi, as shown in Fig. [81 For each value 
of b, beyond some threshold value of Wi, both quantities attain constant values. As a 
consequence, beyond this threshold value, the product Xn^^l^^^* increases linearly with Wi, 
as can be seen clearly in Fig. [TJ enabling a straightforward mapping between Xkxx\x=x* and 
Wi to be made. 

The steep increase in Mxx\x=x* and a*^\x=x* for the Oldroyd-B model, as \nxx\x=x* — ^ 0.5, 
is displayed in Figs. [9] (a) and (b). For the FENE-P model, there is a point of inflection 
at Xkxx\x=x* = 0.5, after which the curves increase much more gradually with increasing 
Xkxx\x=x*- The shapes of the curves for the Oldroyd-B and FENE-P models are strikingly 
reminiscent of the well-known extensional viscosity versus strain rate curves for these mod- 
els, commonly used to disp l ay the unphysical behavior of the Oldroyd-B model 



Bird et al. 



1987b 



Owens and Phillips 



2002| . As is well known, in that case, the onset of the steep 
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Wi Wi 



(a) (b) 

FIG. 8 (Color online) The dependence on Wi of (a) the location x* of the maximum in M^x, and 
(b) the strain rate Kxx\x=x*, for ultra-dilute solutions of Oldroyd-B and FENE-P models. Results 
displayed are solutions of the respective ODEs. 

increase in stress is attributed to the occurrence of a coil-stretch transition, leading to 
an unbounded stress in the Oldroyd-B model, but a bounded stress in the FENE-P model. 
While the former is because of the infinite extensibility of the Hookean spring in the Hookean 
dumbbell model, the latter is because of the existence of a upper bound to the mean stretch- 
ability of the spring in the FENE-P model. We can conjecture, consequently, that in the 
present instance also, polymer molecules undergo a coil-stretch transition in the wake of the 
cylinder, at the location of the stress maximum, giving rise to a stress that increases without 
bound as Wi increases. 

The breakdown of the FEM computations for the Oldroyd-B model is clearly related to 
the steep increase in Mxx\x=x'' and o'xx\x=x*, as \nxx\x=x* — ^ 0.5. (Fig. [9] (b) indicates that 
the stress maximum increases by five orders of magnitude as Wi increases from 0.1 to 4). 
Since the exact solution is known at x = x*, we can calculate the error at any value of Wi 
using, 

error = ^ ^t,^ x 100 (22) 

Error estimates obtained in this manner are displayed in Fig. (a) for the Oldroyd-B 
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(a) (b) 

FIG. 9 (Color online) Coil-stretch transition in the cylinder wake for an ultra-dilute solution. 
Dependence of the polymer stretch (M^x) and stress (o"*^,), at x = x*, on Xkxx\x=x*- The lines are 
ODE solutions and the symbols are FEM results on mesh M5. 

model. The error remains small (< 1%) at low Wi, but increases sharply to approximately 
8% as Wi > The value of Wi*, the Weissenberg number up to which the error is less 

than 1%, depends on the degree of mesh refinement, and a clear improvement in Wi* can 
be observed with increased mesh refinement. However, to obtain mesh converged results for 
Wi > 0.7 (with error < 1%), an approximately 100 fold increase in mesh density and hence, 
approximately ~ 100 fold increase in computational time is required. 

In the case of the FENE-P model, the error in Mxx\x=x* is small even at relatively large 
values of Wi since the chains are close to their fully extended length, and it is difficult to 
see a clear pattern in the change in error with mesh refinement, unlike in the case of the 
Oldroyd-B model above. However, the error in the Myy\x=x* component reveals a more 
systematic behavior, as displayed in Fig. [TD] (b), with a decrease in error with increasing 
mesh refinement. 

The infeasibility of carrying out FEM computations for the confined fiow around a cylin- 
der of an ultra-dilute Oldroyd-B fiuid, for Wi > 1, is revealed in Fig. [TTl where the expo- 
nential increase in the number of elements required to attain an error less than 1%, with 
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FIG. 10 (Color online) Percentage error in FEM computations of (a) M2;a;|a;=a;* for an ultra-dilute 
Oldroyd-B solution, and (b) M.yy\x=x* for an ultra-dilute FENE-P fluid, as a function of Wi. The 
error is calculated based on the results of the ODE solutions. 

increasing Wi, can be clearly observed. For the FENE-P model, the rate of increase in the 
number of elements required for an error less than 1% is significantly lower than that for the 
Oldroyd-B fluid. The Mesh 4 curve in Fig. (TD] (b) even seems to suggest that there might be 
a degree of mesh refinement beyond which, for the FENE-P model, one can compute at any 
Wi with an error less than 1%. However, this trend is not easily discernible in Fig. [11], and 
addtional mesh refinement may be required before a firm conclusion can be drawn. Further, 
a change of variable to the matrix logarithm of the conformation tensor, may lead to mesh 
converged results at significantly higher vales of Wi. 

The reason polymer molecules undergo a coil-stretch transition as they travel down the 
centerline in the cylinder wake is because of the extended period of time they spend in the 
neighborhood of the stagnation point, which leads to a significant accumulation of strain. 
The Hencky strain e at any instant t, calculated from the expression 6 — dt' n^x ) ? is 
displayed in Fig. [121 A very large Hencky strain of roughly 8 units is built up by the time the 
molecules approach x* . For a material element of unit length at t = 0, this corresponds to a 
ratio of final to initial length of roughly 3000. The behavior of individual molecules as they 
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FIG. 11 (Color online) The number of elements as a function of the Weissenberg number Wi* ^ 
up to which the error remains less than 1%, for the Oldroyd-B and FENE-P models. The inset 
displays the same data as a log-log plot. 

are subjected to this degree of straining can be obtained, for an ultra-dilute solution, from 
the Brownian dynamics simulations carried our here since one can calculate the trajectories 
of dumbbells as they are convected by the flow field down the centerline, subjected to the 
local strain rate. 

In this context, it is instructive to calculate the size of individual dumbbells relative 
to a macroscopic feature, such as the length of an element in the finite element mesh. 
In the non-dimensionalization scheme used here, however, since the macroscopic length 
scale is the cylinder radius a, and the microscopic length scale for Hookean dumbbells is 

t unless one ha s estim ates of these length scales. 



(5^)^q/3, direct comparison is difficu 
Here, we use the experimental data of 



McKinley et al. 



1993| . who investigated the flow 



around a confined cylinder (with radius 3.188 x 10 ^m) of a 1.2 million molecular weight 
Polyisobutelene solution, to obtain a typical estimate of these length scales. For these 



molecules, the equilibrium size can be shown to be J{Q^)^^ = 0.0497 /xm. Defining a 
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FIG. 12 (Color online) Hencky strain accumulated by a fluid element as it travels along the 
centerline, from the neighborhood of the stagnation point in the cylinder wake. 

dimensionless length in the flow direction by, 



where, Lm\x=x* is the non-dimensional length of the element at x = x*, the relative length 
of individual molecules in the flow direction can be calculated from the BDS trajectories as 
a function of strain. 

Figures [13] display Q* as a function of e, for an ensemble of 100 dumbbell trajectories, 
at various values of Wi. Nearly all the dumbbells appear to remain close to their initial 
state of extension until approximately 5 strain units, beyond which several of the dumbbells 
undergo rapid extension, which is more pronounced as Wi increases. The rapid extension 
of the dumbbell spring represents the physical unraveling of a polymer molecule from a 
coiled to a stretched state, and the results in Figs. [I3] are inline with the notion that a 
coil-stretch transition occurs as the molecules experience the maximum strain. As expected, 
the molecules relax back to their equilibrium configurations once the strain rate downstream 
of the maximum becomes zero. 

The use of the local element size to achieve non-dimensionalization reveals strikingly that 
the magnitude of some of the molecules is large enough to span several elements. Kinetic 



q: 




(23) 
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m \x=x* 
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FIG. 13 (Color online) The length of individual polymer molecules relative to the size of an 
element of the M5 mesh, at x = x*, at various values of Wi, for an ultra-dilute Oldroyd-B fluid. 



The experimental data of 



McKinlev et al. 



19931 ] is used to obtain an estimate of the equilibrium 



size of the molecules and the cylinder radius. 



theory models, such as the Hookean dumbbell model, are typically built on the assumption 
of homogenous fields, with negligible variation on the length scale of individual molecules. 
The data in Figs. [T3] suggests that the extensive use of the unphysical Oldroyd-B model 
in complex flow simulations is questionable, and highlights the need to derive more refined 
models that are valid in non-homogeneous fields. 



B. Dilute solutions 



All the results reported so far have been for ultra-dilute models, where the existence of 
exact solutions has enabled us to obtain a variety of insights into the origin of difficulties 
encountered with FEM computations. There is already an extensive literature on the nu- 
merical computation of the flow around a confined cylinder of various dilute solution models. 
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FIG. 14 (Color online) Dependence of \kxx\x=x* on Wi for dilute and ultra-dilute solutions. The 
dilute solution curve first approaches the critical value at a value of Wi at which difficulties with 
FEM computations usually arise. 

and there are no new numerical techniques introduced in this paper for us to be able to re- 
port an improvement in the maximum attainable mesh converged Weissenberg number. Our 
interest here, instead, is to examine if any of the insight that has been gained for ultra-dilute 
solutions can be used to understand the observed behavior of dilute solutions. 

The coupling of the velocity and the conformation tensor (and stress) fields, makes it 
impossible to obtain exact solutions for these variables. As a result, it is not possible to 
obtain error estimates as in the case of ultra-dilute solutions, or to calculate the trajectories 
of individual dumbbell molecules convected along the centerline by the flow. Nevertheless, 
we can exploit the key insight of the previous section because the value of the maximum in 
the component in the cylinder wake, for the Oldroyd-B model, is still given by eqn. (1211) . 
even though we do not have a pre- determined velocity field Vx{x). 

In the case of an ultra-dilute solution, it was relatively straightforward to find the location 
X* of the stress maximum in the cylinder wake for each value of Wi, and to calculate Xnxx\x=x* 
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FIG. 15 (Color online) Coil-stretch transition in the cylinder wake for a dilute solution. Variation of 
the polymer stretch M^^; with Xkxx\x=x* for Oldroyd-B and FENE-P models. The M4 mesh is used 
for the Oldroyd-B model and the M3 mesh is used for the FENE-P model in FEM computations. 
The viscosity ratio (3 = 0.59. 

from the known velocity field. For a dilute solution, the velocity field changes with each 
change in Wi. All the same, it is still possible to find \nxx\x=x''i at various values of Wi, by 
integrating the full system of equations for the velocity and conformation tensor fields using 
the FEM formulation. A typical velocity profile for a dilute solution at Wi = 0.6 is shown 
in Fig. [3] as the dot-dashed line. (At this value of Wi, it can be seen that the velocity profile 
is not significantly different from that for an ultra-dilute solution.) 

Figure UM indicates that Xk,xx\x=x* ioi a dilute solution approaches the critical value of 0.5 
in a manner similar to that of an ultra-dilute solution, albeit at a slightly more rapid rate with 
increasing Wi. The curve is not as smooth as the ultra-dilute case because of the growing 
error in FEM computations as Wi reaches the upper limit of computable values. Indeed, the 
typical limit of Wi = 0.6 where most computations in the literature first encounter problems, 
appears to be the value at which Xkxx\x=x* first comes close to 0.5. 



31 



As may be anticipated from eqn. (12T|) . M^x and a*^\x=x* will increase steeply for the 
Oldroyd-B model as Xkxx\x=x* —>■ 0.5. This can be seen very clearly from Fig. [151 where 
Mxx appears to become unbounded in this limit. For the FENE-P fluid on the other hand, 
the curves for M^^^^ exhibit a point of inflection at Xk^^I^^^^* = 0.5, before leveling off to the 
fully stretched value corresponding to the respective value of b. 

The similarity of the shapes of the curves in Fig. [15] to the curves in Fig. [9l and to the 
well-known extensional viscosity versus strain rate curves for the Oldroyd-B and FENE-P 
models suggests that even for a dilute solution, there occurs a coil-stretch transition at the 
location of the stress maximum in the cylinder wake, and this coil-stretch transition is the 
source of problems encountered with FEM computations for the flow of an Oldroyd-B fluid 
around a confined cylinder. It would be of great interest to examine if a similar coil-stretch 
transition is the source of computational difficulties encountered in the numerical simulation 
of other benchmark complex flows of Oldroyd-B fluids. 

V. CONCLUSIONS 

The flow around a cylinder confined between parallel plates, of ultra-dilute and dilute 
polymer solutions, modeled by the Oldroyd-B and FENE-P constitutive equations, has been 
considered with a view to understand the origin of computational difficulties encountered in 
numerical simulations. 

FEM computations of ultra-dilute Oldroyd-B solutions are shown to breakdown at 
Wi = 0{1), as has been observed previously for dilute solutions (see Fig. [4]). Two dif- 
ferent numerical means of obtaining an exact solution along the centerline in the cylinder 
wake, for both the Oldroyd-B and the FENE-P models, have been developed to enable a 
careful examination of the causes of the breakdown. The exact solution techniques, namely 
solving a system of ODEs and carrying out Brownian dynamics simulations, are useful to 
evaluate the value of Wi up to which the FEM computations are accurate (see Figs. [5]), and 
to estimate the error in the FEM results (see Figs. [TOj) . 

An analysis of the structure of the Oldroyd-B equation shows that the maximum in the 
Mxx component of the conformation tensor along the centerline in the cylinder wake, be- 
comes unbounded if the non-dimensional strain rate at the location of the stress maximum, 
\nxx\x=x*, approaches a critical value of 0.5 (see eqn Numerical solution of the gov- 
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erning equations for an ultra-dilute solution reveals that \kxx\x=x* 0.5 as a power-law in 
Wi, for Wi > 0{1) (see Fig. [6] (a) and (b)). On the other hand, analysis of the FENE-P 
model reveals that the maximum in M^.^. remains bounded for all values of the local non- 
dimensional strain rate. In contrast to the Oldroyd-B model, numerical results for \Kxx\x=x* 
show that it does not have an asymptotic value, but instead increases linearly with Wi 
beyond a threshold value of the Weissenberg number (see Fig. [7]). 

As \nxx\x=x* ^ 0.5, both the maximum in Vixx and the maximum in the stress (o"*^) 
increase without bound for the ultra-dilute Oldroyd-B model. In comparison, for an ultra- 
dilute FENE-P model, these variables increase relatively rapidly as \nxx\x=x* approaches 
0.5, but level off and remain bounded for higher values of (see Figs. IH]). The shape 

of the curves are strongly suggestive of the occurrence of a coil-stretch transition in the 
cylinder wake. 

The steep increase in yixx and a*^ in the vicinity oiWi = 0{1) necessitates the use of 
increasingly refined meshes for increasing values of Wi. The number of elements required to 
maintain the error in Mxx less than 1% is shown to increase exponentially with increasing 
Wi for the Oldroyd-B model, making it practically infeasible to obtain solutions at Wi > 1. 
In the case of the FENE-P model, the current simulation data is inadequate to draw firm 
conclusions (see Figs. [TTj) . 

A material element of an ultra-dilute solution is shown to accumulate nearly 8 units 
of Hencky strain as it travels downstream from the stagnation point in the wake of the 
cylinder, due to the extended time it spends in the vicinity of the stagnation point (see 
Fig. fT2l) . This strain leads to dumbbells undergoing a large extension in the flow direction, 
with their magnitude large enough to span several elements in the local finite element mesh 
see Fig. [13]). 

The analysis of the nature of the maximum in Mxx in the cylinder wake, which suggests 
that the maximum becomes unbounded if Xkxx\x=x* approaches the critical value of 0.5, 
is valid for both ultra-dilute and dilute Oldroyd-B fluids. FEM computations of the fully 
coupled governing equations for a dilute Oldroyd-B fluid have been carried out to show 
that, just as in the case of an ultra-dilute solution, Xkxx\x=x* tends to 0.5 with increasing 
Wi (see Fig. [T^ . The approach, however, is more rapid than the ultra-dilute case, with 
Xkxx\x=x* becoming nearly equal to 0.5 at Wi ~ 0.6, the value of the Weissenberg number 
where computational difficulties have been reported in the literature to be first typically 
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encountered. 

The approach of \kxx\x=x* to the critical value is shown to be accompanied by an un- 
bounded increase in the maximum values of M^.^ and o"*^ in the cylinder wake for a di- 
lute Oldroyd-B fluid. FEM computations of the coupled governing equations for a FENE- 
P fluid on the other hand show that these variables increase relatively rapidly close to 
\tixx\x=x* = 0.5, but level off and remain bounded at higher values (see Figs. [15]). The 
similarity of the curves with observations for ultra-dilute solutions, is strong evidence for a 
coil-stretch transition also occurring in dilute solutions, in the wake of the cylinder at the 
location of the stress maximum. 

Several issues that must be addressed in the future can be tackled fruitfully with the 
framework developed here. For instance, the nature and structure of stress boundary layers 
in the vicinity of the cylinder can be examined for ultra-dilute solutions along lines similar 
to the analysis here. Further, the existence of a coil-stretch transition suggests that a model 
with conformation dependent drag might reveal the existence of coil-stretch hysteresis in the 
cylinder wake. 
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